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The  power  management  strategy  (PMS)  plays  an  important  role  in  the  optimum  design  and  efficient 
utilization  of  hybrid  energy  systems.  The  power  available  from  hybrid  systems  and  the  overall  lifetime  of 
system  components  are  highly  affected  by  PMS.  This  paper  presents  a  novel  method  for  the  determina¬ 
tion  of  the  optimum  PMS  of  hybrid  energy  systems  including  various  generators  and  storage  units.  The 
PMS  optimization  is  integrated  with  the  sizing  procedure  of  the  hybrid  system.  The  method  is  tested  on  a 
system  with  several  widely  used  generators  in  off-grid  systems,  including  wind  turbines,  PV  panels,  fuel 
cells,  electrolyzers,  hydrogen  tanks,  batteries,  and  diesel  generators.  The  aim  of  the  optimization  problem 
is  to  simultaneously  minimize  the  overall  cost  of  the  system,  unmet  load,  and  fuel  emission  considering 
the  uncertainties  associated  with  renewable  energy  sources  (RES).  These  uncertainties  are  modeled  by 
using  various  possible  scenarios  for  wind  speed  and  solar  irradiation  based  on  Weibull  and  Beta  proba¬ 
bility  distribution  functions  (PDF),  respectively.  The  differential  evolution  algorithm  (DEA)  accompanied 
with  fuzzy  technique  is  used  to  handle  the  mixed-integer  nonlinear  multi-objective  optimization  prob¬ 
lem.  The  optimum  solution,  including  design  parameters  of  system  components  and  the  monthly  PMS 
parameters  adapting  climatic  changes  during  a  year,  are  obtained.  Considering  operating  limitations  of 
system  devices,  the  parameters  characterize  the  priority  and  share  of  each  storage  component  for  serving 
the  deficit  energy  or  storing  surplus  energy  both  resulted  from  the  mismatch  of  power  between  load  and 
generation.  In  order  to  have  efficient  power  exploitation  from  RES,  the  optimum  monthly  tilt  angles  of 
PV  panels  and  the  optimum  tower  height  for  wind  turbines  are  calculated.  Numerical  results  are  com¬ 
pared  with  the  results  of  optimal  sizing  assuming  pre-defined  PMS  without  using  the  proposed  power 
management  optimization  method.  The  comparative  results  present  the  efficacy  and  capability  of  the 
proposed  method  for  hybrid  energy  systems. 

©  2011  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Distributed  generation  (DG)  has  been  recently  nominated  as 
one  of  the  solutions  to  reliable,  less  costly,  and  more  efficient 
energy  supply  systems.  Specifically,  small  DG  systems  with  power 
level  ranging  from  1  kW  to  10  MW  [1],  located  near  the  loads, 
are  extensively  utilized  both  in  grid-connected  and  stand-alone 
modes.  Among  available  DGs,  renewable-based  systems  (RES) 
such  as  photovoltaics  and  wind  turbines  have  attained  the  most 
popularity  due  to  ever  increasing  concerns  about  depletion  of 
fossil  fuels  and  global  warming.  They  have  also  been  getting  more 
cost-effective  during  recent  years. 

However,  a  significant  drawback  associated  with  solar  and  wind 
energy  systems  is  their  intermittent  and  unpredictable  behavior 
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due  to  their  direct  dependence  on  climatic  conditions.  In  addition, 
the  variations  of  the  wind  and  solar  energy  may  not  match  with  load 
changes  [2].  The  reliability  of  the  system  is  important  to  both  plan¬ 
ning  and  utilization  stages.  Designing  energy  systems  including 
solar  and  wind  energies  together,  to  some  extent,  reduces  the  depth 
of  the  problem.  Since,  the  advantage  of  one  source  can  overcome 
the  disadvantage  of  the  other  one  and  vice  versa  [3].  In  addition, 
taking  into  account  the  intermittency  and  uncertainty  associated 
with  solar  and  wind  energy  sources,  improves  the  adaptation  of 
design  results  with  practical  and  realistic  conditions.  Another  con¬ 
ventional  solution  is  appropriate  incorporation  of  energy  storage 
devices  such  as  batteries  and  hydrogen  storage  systems  compris¬ 
ing  electrolyzers,  hydrogen  tanks  and  fuel  cells  into  the  system 
[3].  Hence,  hybrid  energy  systems  are  becoming  more  attractive 
to  power  engineers. 

One  of  the  most  prominent  issues  regarding  hybrid  energy 
systems  is  to  determine  their  optimum  design  and  operation 
modes  taking  account  of  regional  conditions  and  load  demand 
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characteristics.  Therefore,  many  efforts  have  been  made  for  design¬ 
ing  and  planning  hybrid  DG  systems  [1-18].  These  works  can  be 
reviewed  from  various  points  of  view  such  as  system  modeling,  the 
applied  data  in  the  solution  algorithm,  objective  functions  to  be 
optimized  and  the  mathematical  tool  used  to  handle  the  optimiza¬ 
tion  problem.  In  [2],  the  overall  cost  of  the  introduced  hybrid  system 
throughout  the  estimated  lifespan  of  the  installation  has  been 
minimized.  The  reliability  of  serving  energy  and  proper  functioning 
of  components  has  been  involved  in  the  applied  algorithm  as  a  con¬ 
stant  constraint  not  to  be  violated.  Solar  radiation  and  wind  speed 
are  also  assumed  deterministic.  Particle  Swarm  Optimization  (PSO) 
algorithm  has  been  applied  to  handle  the  single-objective  problem. 
In  [3],  the  cost  analysis  of  a  hybrid  wind/PV/fuel  cell  system  has 
been  focused  on  through  a  residential  consumer  case  study. 

In  [4-10],  the  cost  has  been  minimized  by  selecting  the  corre¬ 
sponding  system  configuration  and/or  operation  strategy.  The  load 
demand  has  been  supposed  to  be  completely  or  partially  fulfilled 
with  a  fixed  user-defined  value  for  unmet  load.  In  this  case,  the 
solution  method  has  no  choice  to  reach  a  compromise  between 
cost  and  unmet  load.  In  [9,10],  pollutant  emissions  have  been  rep¬ 
resented  in  the  form  of  cost  and  limited  by  being  economically 
involved  in  the  cost  function,  and  thus  the  results  depend  on  the 
cost  coefficient  assigned  to  the  emissions.  More  realistic  optimal 
configuration  may  be  achieved  by  involving  pollutant  emissions 
measurement  in  their  relative  unit  (kg  per  liters  of  fuel)  instead  of 
converting  it  to  the  cost.  From  the  aspect  of  the  system  modeling 
accuracy,  Ref.  [5],  has  taken  a  fixed  PV  panels  tilt  angle  into  account 
using  Hay  and  Devis  (1980)  and  Orgin  and  Holland  (1977)  meth¬ 
ods,  whereas  [6]  has  imported  the  angle  value  as  an  optimization 
variable  in  the  solution  algorithm.  In  the  latter  work,  commercially 
available  types  for  each  device  and  their  costs  have  been  taken 
into  account,  and  optimization  results  representing  size  and  type 
of  each  device  have  been  demonstrated.  The  load  demand  has  been 
assumed  to  be  entirely  fulfilled.  The  maximum  power  point  tracker 
(MPPT)  has  been  included  in  both  PV  and  wind  turbines  systems. 
The  presented  approach  in  [7]  has  considered  the  influence  of  the 
uncertainty  in  solar  irradiation  on  the  sizing  process  of  an  off-shore 
PV  power  system  using  three  probabilistic  models  such  as  Markov 
Chain.  A  size  optimization  procedure  has  been  demonstrated  in 
[8],  where  the  results  have  been  compared  by  using  the  optimiza¬ 
tion  software  called  HOMER  [19].  It  has  been  stated  that  because  of 
the  nonlinearity  and  complexity  of  hybrid  systems,  the  application 
of  evolutionary  algorithms  like  genetic  algorithm  generates  better 
results  than  the  application  of  classical  optimization  methods  [8]. 

On  the  case  of  PMS  [5],  as  one  of  the  earliest  works,  has 
introduced  two  parameters  named  SDM  and  SAR,  representing 
the  minimum  and  maximum  State  of  Charge  (SOC)  of  batter¬ 
ies,  respectively.  The  optimum  values  of  the  SOCs  have  been 
calculated  within  the  context  of  the  size  optimization  problem. 
Although  [6]  has  applied  an  accurate  model  for  PV  and  wind 
turbines  systems  representing  installation  parameters,  the  oper¬ 
ation  control  has  been  fixedly  defined  by  the  user  and  did  not 
conform  with  the  system  conditions.  In  [8],  the  effect  of  the  min¬ 
imum  and  maximum  limit  of  SOC  on  the  system  operation  and 
cost  has  been  evaluated.  The  algorithm  suffers  a  computation  bur¬ 
den  since  a  sub-algorithm  for  finding  optimum  values  of  control 
parameters  has  been  executed  for  any  single  control  vector  gener¬ 
ated  in  the  main  optimization  algorithm  as  a  nested  optimization 
loop. 

In  [9],  a  grid-connected  hybrid  energy  system,  capable  of 
mutual  exchange  of  power  with  the  grid,  comprising  wind  turbine, 
micro-turbine,  diesel  generator,  photovoltaic  array,  fuel  cell  and 
battery  storage  has  been  analyzed  and  on-line  optimal  power 
management  has  been  attained.  The  problem  has  been  treated  as 
a  single  objective  problem  considering  all  objectives  such  as  fuel 
emissions  in  terms  of  cost.  The  fuel  emission  has  been  factorized 


into  three  main  gas  components  and  values  for  each  gas  has  been 
separately  presented. 

The  cost  minimization  of  a  hybrid  energy  system  including 
hydrogen  storage  has  been  presented  in  [10].  Using  the  classic 
linear  programming  optimization  tool,  a  comparison  between  an 
optimized  dispatch  strategy  and  a  fixed  one  has  been  performed  so 
that  the  latter  strategy  can  be  reformed  using  the  values  obtained 
for  the  dispatch  variables  of  the  former  one.  Ref.  [11]  has  incorpo¬ 
rated  the  effect  of  the  number  of  battery  charge/discharge  cycles 
relative  to  their  depth  of  discharge  (DOD)  on  their  lifespan,  oper¬ 
ation  strategy  and  consequently  the  system  cost.  A  similar  work 
has  been  addressed  in  [12],  which  has  modeled  the  uncertainty 
involved  in  operating  and  design  characteristics  of  the  system  tak¬ 
ing  the  advantage  of  Stochastic  Annealing  optimization  algorithm. 
This  work  has  referenced  the  classification  of  uncertainties  pre¬ 
sented  in  [13]  as  model-inherent,  system  inherent  and  external 
uncertainties  existing  in  the  system. 

In  [14],  a  method  has  been  introduced,  in  which  the  optimized 
configuration  and  operation  of  the  system  has  been  achieved 
based  on  twelve  parameters  defined  relative  to  operational  limits 
of  different  storage  devices.  This  work  has  been  rearranged  in 
[15]  considering  a  more  complex  hybrid  system  and  a  couple  of 
objective  functions  including  cost,  emission  and  unmet  load. 

As  seen  in  [14],  the  power  available  from  such  systems  and 
the  overall  lifetime  of  system  components  are  highly  affected  by 
the  applied  power  management  strategy;  hence,  various  strate¬ 
gies  result  in  various  designs  for  the  system  aiming  to  meet  load 
requirements  within  the  estimated  lifespan  of  the  system. 

In  this  work,  a  general  method  for  optimal  PMS  of  autonomous 
systems  including  various  generators  and  storage  media,  com¬ 
bined  with  optimal  design  of  system  components  is  proposed. 
To  do  so,  appropriate  model  of  system  components  for  the  plan¬ 
ning  problem  is  reviewed.  To  solve  the  introduced  problem,  the 
differential  evolution  algorithm  (DEA),  which  is  nominated  to  be 
capable  of  solving  various  nonlinear  optimization  problems,  is  uti¬ 
lized.  The  algorithm  is  accompanied  with  the  fuzzy  technique  to 
handle  the  multi-objective  optimization  problem  in  a  time  effi¬ 
cient  manner.  In  this  case,  there  is  no  need  to  find  appropriate 
coefficients  as  penalty  factors  of  constraints  violations  and  as 
objective  function  weights.  A  novel  method,  involved  in  the  size 
optimization  problem,  is  proposed  to  obtain  values  for  the  param¬ 
eters  of  PMS  by  which  the  optimum  performance  and  minimum 
cost,  emission  and  unmet  load  are  achieved.  Considering  operat¬ 
ing  limitations  of  system  devices,  these  parameters  characterize 
the  priority  and  share  of  each  storage  component  for  serving 
the  deficit  energy  or  storing  surplus  energy  both  resulted  from 
the  mismatch  of  power  between  load  and  RES.  Optimal  values 
for  design  parameters  and  PMS  parameters  are  simultaneously 
attained. 

The  RES  uncertainties  are  applied  to  the  optimization  procedure 
by  scenario  generation  based  on  Weibull  [20]  and  Beta  distribu¬ 
tion  functions  for  wind  speed  and  solar  irradiation,  respectively 
[21].  Numerical  results,  including  type  and  number  of  each  com¬ 
ponent,  monthly  values  for  PV  panels  tilt  angle,  the  height  of  wind 
turbine  towers  along  with  the  PMS  parameters  ensure  the  capabil¬ 
ity  of  the  proposed  method  to  achieve  the  aim  of  optimization.  In 
order  to  focus  on  the  significance  of  PMS  and  show  the  influence  of 
strategy  variations  on  the  performance  and  objectives  of  the  sys¬ 
tem,  a  comparison  is  presented,  too.  In  addition,  determining  the 
PMS  parameters  in  the  optimization  procedure  is  confirmed  as  an 
efficacious  concept. 

In  the  next  section,  the  description  of  the  system  under  study 
and  the  model  of  components  are  presented.  Section  3  illustrates 
the  problem  statement  including  objective  functions,  constraints 
and  the  proposed  method  for  the  power  management.  In  Section 
4,  DEA  and  fuzzy  multi-objective  technique  are  explained.  Finally, 
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the  results  and  conclusions  are  represented  in  two  subsequent 
sections. 

2.  Model  of  system  components 

In  this  work,  a  stand-alone  hybrid  energy  system  comprising 
wind  turbines,  PV  panels,  fuel  cells,  electrolyzers,  hydrogen  tanks, 
batteries  and  diesel  generators  is  studied.  A  simplified  diagram  of 
the  system  is  depicted  in  Fig.  1.  The  available  power  from  renew¬ 
able  sources  is  directly  delivered  to  the  load  in  order  to  provide  the 
load  demand.  The  excess  or  deficient  power  from  RES  is  saved  in  or 
supplied  from  other  equipments  in  the  system  based  on  the  pro¬ 
posed  dispatch  strategy  explained  in  Section  3.2.  The  model  of  each 
system  component  is  briefly  described  in  the  following  subsections. 

2.2.  Wind  turbine 

The  output  power  of  the  wind  turbine,  Pw,  is  calculated  on  the 
basis  of  its  power  curve  given  by  the  manufacturer.  The  effect  of 
the  wind  turbine  installation  height  and  the  roughness  of  the  land 
surface  of  the  towers  on  the  output  power  is  also  considered  [2,6]. 

2.2.  PV  panel 


The  solar  elevation  angle  (or),  which  is  the  angle  between  the  direc¬ 
tion  of  the  sun  and  the  horizon  is  then  estimated,  as  follows: 

o  =  sin* 1  [sin  (p  sin  8  +  cos  cp  cos  <5  cos  (LMST  -  A.)]  (2) 

where  LMST  stands  for  Local  Mean  Sidereal  Time  and  cp  is  the  geog¬ 
raphy  of  the  longitude  [17].  The  diffuse  and  reflected  radiations  are 
neglected.  Using  Fig.  2,  the  incident  radiation  on  the  tilted  surface 
(Gj)  can  be  expressed  in  terms  of  horizontal  component  of  solar 
irradiation  ( Gh ),  as  follows: 


-c 

o 

II 

o 

(3) 

sin  a 

Gp  =  Gj  sin(<7  +  /3) 

(4) 

where  Gp  is  the  effective  component  of  the  solar  radiation  perpen¬ 
dicular  to  the  panel.  Based  on  the  calculated  Gp,  the  power  available 
from  PV  panels  at  time  step  t  considering  the  temperature  effect  is 
determined  using  the  following  equations: 

NCOT  -  20 

Tdt )  =  TA(t )  +  8QQ  Cp(t,  fi)  (5) 

fec(0  =  Usc,stc  +  I<i(Tc(t)  -  25°)]  (6) 

Voc(t)  =  V0c,stc  -Kv-Tc(t)  (7) 

PM(t)  =  Ns-Np.Voc-ISdt)-FF(t)  (8) 

where  NCOT  is  the  nominal  cell  operating  temperature  (°C),  TA(t) 
is  the  ambient  temperature  at  t,  Tc(t)  is  the  cell  temperature  at  t, 
V0c,stc  and  Isc,stc  are  the  module  open-circuit  voltage  and  short- 
circuit  current  under  STC,  with  K\  and  I<y  as  their  corresponding 
temperature  coefficients;  Pm(0  is  the  power  of  array  with  NP  mod¬ 
ules  in  parallel  and  Ns  modules  in  series  and  FF{t)  is  the  fill  factor 
[6]. 

2.3.  Fuel  cell 


To  achieve  the  maximum  output  power  from  a  PV  panel,  it  is 
essential  to  consider  the  influence  of  the  solar  position  and  panel  tilt 
angle.  In  [  1 6],  it  is  shown  that  such  effects  can  change  the  total  avail¬ 
able  power  up  to  9.94%  of  the  total  maximum  power.  The  panels 
are  assumed  to  have  rotation  on  one  axis.  The  following  equations 
lead  to  the  incident  irradiation  on  a  tilted  panel.  Eq.  (1 )  is  employed 
to  calculate  the  longitude  at  the  equator  (<$),  with  respect  to  solar 
ecliptic  longitude  and  latitude  (A,  0)  and  earth’s  axis  inclination  to 
the  plane  of  its  orbit  (/  =  23.439  °). 


The  formulations,  which  represent  the  balance  between  the 
hydrogen  as  input  and  the  electric  power  as  output  in  steady  state 
condition,  is  the  required  fuel  cell  model  used  in  the  design  pro¬ 
cedure.  However,  for  the  sake  of  ensuring  maximal  efficiency  of 
the  system  and  preventing  degradation  of  fuel  cells  due  to  abun¬ 
dance  of  start-up  and  shut-down  cycles,  a  hysteresis  band  is  also 
embedded  into  the  model,  as  in  [12]. 

The  fuel  cell  hydrogen  consumption,  ConsFC  (kg/h),  relative  to 
its  output  electric  power  is  computed  by  the  following  equation 
[15]: 


Cons  Fc  = 


Bfc  •  Pn_fc  +  Ape  ■  Pfc 

Bfc  ■  Pn_fc  +  AFc  •  Pfc  +  Fef 


Pmax  _e/ 


forP/PN_FC 

—  Pmax  _ef 

otherwise 


<5  =  sin  j(sin  0  cos  I  +  cos  0  sin  X  sin  I) 


where  PFc  is  the  output  power  of  fuel  cell  in  kW,  PNFc  is  the  nomi- 

(1 )  nal  output  power,  AFC  and  BFC  (kg/kWh)  are  the  coefficients  of  the 
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consumption  curve.  Pmax_e/>  in  terms  of  percentage  of  Pn_fc,  is  the 
output  power,  at  which  the  efficiency  of  the  fuel  cell  is  maxi¬ 
mum  and  Fef  is  a  factor  to  consider  the  high  consumption  above 
Pma xjef-  The  hydrogen  consumption  parameters  for  all  fuel  cells  are 
assumed  to  be  AFG  =  0.05  and  BFC  =  0.004  kg/ kWh,  Pmax  ef  =  0.2  and 
Per  1  [15].  By  using  these  values,  the  fuel  cell  efficiency  is  46%  at 
the  maximum  and  31%  at  rated  powers. 

2.4.  Electrolyzer  and  hydrogen  tank 

The  input  electrical  energy  dependence  on  the  hydrogen  mass 
flow  is  modeled,  as  follows  [15]: 

ConsE  =  BE  ■  Qn _e  +  AE  ■  Q  (10) 


where  AG  and  BG  (1/kW)  are  fuel  consumption  curve  coefficients 
provided  by  the  manufacturer,  PG  (kW)  is  the  output  power  and 
PN  G  (1<W)  is  the  nominal  output  power  of  the  diesel  generator. 

The  values  assigned  to  AG  and  BG  are  0.246  and  0.08145  1/kWh, 
respectively,  for  all  diesel  generators,  which  have  been  used  in  [  1 5], 
too. 

3.  Problem  statement 

3.1.  Objective  functions  and  constraints 

The  objective  functions,  considered  in  the  optimization  prob¬ 
lem,  are  as  follows: 


where  Qn_e  is  the  nominal  hydrogen  mass  flow  (kg/h),  Q.  is  the 
hydrogen  mass  flow  (kg/h),  AF  and  BE  are  the  coefficients  of  the 
consumption  curve  in  kW/kg/h.  The  parameters  of  the  electrical 
energy  consumption  for  all  electrolyzers  are  assumed  to  be  the 
same  as  parameters  given  in  [15],  i.e.,  AF  =  40  and  BF  =  20  both  in 
kW/kg/h. 

The  maximum  hydrogen  tank  capacity  is  assumed  to  be  equal  to 
10  kg  and  the  total  capacity  of  the  hydrogen  storage  is  determined 
by  the  design  program  as  the  number  of  hydrogen  tanks  (Ntan/J. 
The  minimum  allowed  level  of  the  hydrogen  (in  tanks)  is  also  a 
parameter  determined  by  the  proposed  PMS. 


2.5.  Battery 


The  battery  bank  with  the  total  nominal  capacity  of  Cn  (Ah),  can 
serve  energy  to  the  load  until  the  maximum  depth  of  discharge 
(DOD)  or  SOCmin  is  reached.  SOCmin  of  battery  bank  along  with 
SOCmax  are  taken  into  account  as  two  of  the  parameters  of  PMS, 
which  should  be  optimized  by  the  developed  program  and  should 
not  exceed  the  values  introduced  by  the  manufacturer. 

The  SOC  of  battery  bank,  in  each  simulation  time  step,  is  calcu¬ 
lated  by  applying  the  following  equation: 


SOC(t)  =  SOC(t-l)  +  nB^llOO  (11) 

where  SOC(t)  is  the  batteries’  SOC  at  time  step  t,  nB  is  the  battery 
round-trip  efficiency  and  PB(t)  is  the  power  charged  in  or  dis¬ 
charged  from  battery  bank  at  time  step  t  [6].  nB  is  approximately 
equal  to  80%  in  charging  and  100%  in  discharging  modes  [6].  PB(t ) 
is  determined  due  to  the  mismatch  of  power  between  load  and 
RES  and  the  share  associated  with  batteries  to  supply  the  defi¬ 
cient  power  or  to  save  the  excess  power.  The  mentioned  share  is  a 
parameter  of  PMS  and  optimized  in  the  proposed  algorithm. 

The  number  of  series  connected  batteries  depends  on  the  DC  bus 
voltage  ( VBUs )  and  the  nominal  voltage  of  each  battery,  i.e.: 


,S  _  VBUS 

B~  VB 


(12) 


The  total  number  of  batteries  in  battery  bank  is  equal  to  n|  multi¬ 
plied  by  (the  number  of  batteries  in  parallel),  which  should  be 
determined  in  the  optimization  algorithm. 

In  order  to  model  the  effect  of  the  number  of  charge/discharge 
cycles  on  the  lifetime  of  batteries,  the  model  used  in  HOMER  soft¬ 
ware  [22]  is  used.  This  model  is  based  on  the  battery  “cycles  to 
failure”  curve,  which  is  more  explained  in  detail  in  [1 1  ]. 


2.6.  Diesel  generator 

The  fuel  consumption  of  the  diesel  generator,  ConsG,  in  terms  of 
its  output  power  is  written,  as  follows: 


ConsG  =  Bg  •  P]v_g  +Ag-Pg 


(13) 


•  The  overall  cost  of  the  system  discounted  to  the  year  of 
installation  (NPC)  including  the  investment  cost,  operation  and 
maintenance  cost  during  the  lifespan  of  the  system,  replacement 
cost,  and  fuel  cost  of  diesel  generators.  All  of  the  cost  terms  asso¬ 
ciated  with  the  system  components  are  based  on  the  data  in  [  1 5]. 

•  The  total  fuel  emissions  produced  by  diesel  generators  during  the 
total  lifespan  of  the  system.  Thanks  to  the  employed  fuzzy  multi¬ 
objective  optimization  technique,  the  model  implemented  for  the 
fuel  emission  is  adequately  accurate  and  there  is  no  requirement 
for  taking  it  into  account  in  terms  of  cost.  As  the  most  significant 
gas  existing  in  the  diesel  generator  exhaust  is  C02,  which  also 
causes  green  house  effects,  the  produced  C02  is  considered  to 
represent  the  fuel  emission. 

•  The  total  Loss  of  Load  Probability  (LPSP),  which  represents  the 
total  unmet  energy.  This  index  is  the  most  commonly  used  index 
in  related  works  [2].  LPSP  is  proportional  to  the  unmet  load  and 
is  written,  as  follows: 


LPSP  = 


LOEE 


(14) 


where  D(/i)  is  the  load  demand  in  kWh  at  time  step  h,  and  LOEE , 
standing  for  Loss  of  Energy  Expected,  is  defined  by  the  following 
equation: 

Eld 

E[LOE(h)]  (15) 

h= l 


E[LOE]  =  ^Q(s)-/(s)  (16) 

sgS 

where  Q(s)  represents  the  amount  of  the  unmet  energy  when  the 
system  experiences  the  state  s  and/(s)  is  the  probability  associated 
with  the  occurrence  of  the  state  s. 

The  control  variables,  which  should  be  used  in  simulations  to 
calculate  the  objective  functions,  as  well  as  their  boundaries,  are 
listed,  as  follows: 


-  A  k  x  2  matrix  [[Nk],  [Tk]]  of  size  parameters  of  the  system  com¬ 
ponents  shown  in  Fig.  1,  where,  Nk  is  the  number  and  Tk  is  the 
type  of  component  k.  Both  are  positive  integers. 

-  A  vector  of  installation  parameters  including  monthly  tilt  angles 
of  PV  panels  (0°  </3/<  <90°,  /<=  1,. .  .,12)  and  wind  turbine  tower 
height  (h>  0), 

-  A  vector  of  PMS  parameters  (described  in  the  coming  section) 
including: 

-  Monthly  charge  shares  associated  with  n  -  1  types  of  storage 
media  (0<CSj<l,  where,  n  is  the  total  types  of  the  storage 
media  existing  in  the  system), 

-  Monthly  discharge  shares  associated  with  m  -  1  types  of  backup 
devices  (0<DSZ  <1,  where,  m  is  the  total  number  of  backup 
devices  capable  of  feeding  power  to  the  load), 
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-  Upper  and  lower  limits  on  the  stored  energy  level 
of  all  the  storage  media  (50Q  min  rated  <SOQ  min  <  0.5, 
0.5  <  SOCi  max  <  S0Q>max  ratecf,  0  <  SOC  <  1 ). 

In  addition  to  the  boundaries  that  control  variables  should  not 
take  values  beyond  them,  there  are  operational  constraints  to  be 
mentioned  as  follows: 


The  value  of  SOCmax  is  obtained  within  the  context  of  the  optimiza¬ 
tion  procedure. 

Case  ( 1  -3 ):  P,(  t)  is  beyond  the  HBR.  Hence,  the  power  to  be  drawn 
from  the  device  does  not  fall  within  the  allowable  range  and  no 
power  can  be  delivered,  i.e.: 

P*(t)  =  0  (21) 


-  The  specified  maximum  and  minimum  SOC  of  batteries  (SOCmin, 
SOCmax  )> 

-  The  rated  current  or  power  of  all  devices, 

-  The  minimum  power  injected  to  the  electrolyzers  and  the  min¬ 
imum  supplied  power  from  fuel  cells  and  diesel  generators,  at 
which  they  can  start-up.  This  constraint  is  considered  to  avoid 
the  multitude  of  the  start-up/shut-down  cycles  or  short  periods 
of  frequent  utilization  of  these  devices,  and  enhance  their  life¬ 
time.  This  operation  limit  has  been  introduced  as  Hysteresis  Band 
Range  (HBR)  in  [12].  Other  imposed  constraints  are  illustrated  in 
the  next  section. 

3.2.  Proposed  method  for  power  management  of  system 

In  general,  the  PMS  is  so  important  that  the  operation,  relia¬ 
bility,  cost  and  lifetime  of  the  system  is  affected  with  even  minor 
alterations  in  the  strategy.  In  the  proposed  dispatch  strategy,  to  effi¬ 
ciently  utilize  the  RES,  the  available  power  from  RES,  i.e.,  PresW,  is 
directly  delivered  to  the  load  in  order  to  provide  the  load  demand, 
P/oad(t).  Other  equipments  are  considered  without  any  predefined 
priorities  and  there  is  not  any  preference  to  some  components.  The 
priorities  are  determined  by  the  optimization  process.  The  excess 
or  deficient  power  of  RES,  i.e.,  Prem(t),  is  saved  in  or  supplied  from 
other  equipments  in  the  system,  i.e.: 


Prem(f)  =  PrES^)~  PloadW  07) 

where  Prem(t)  might  be  a  negative  or  positive  value. 

For  Prem(t)  >  0:  In  this  case,  PresW  is  more  than  the  load  demand 
at  time  step  t,  and  the  remaining  power  should  be  delivered  to  the 
storage  devices  namely  electrolyzers  and  battery  bank.  Therefore, 
based  on  the  proposed  method,  the  optimization  procedure  should 
assign  Charge  Shares  (CS)  for  battery  bank  and  electrolyzers  so  that 
Prem(t)  is  optimally  shared  among  them  considering  constraints. 
The  power,  to  be  stored  in  storage  device  /,  i.e.,  P,(t),  is  written,  as 
follows: 

Pft)  =  CSfk)  ■  Prem(t)  i=l,...,n  (18) 

where  n  is  the  total  number  of  storage  media  and  CSfk)  is  the  CS 

of  storage  device  i  in  month  k  (/<  =  1 . 12).  It  should  be  mentioned 

that  all  CSs  fall  within  the  range  (0,1 ). 

Now,  the  status  of  all  storage  devices  is  checked  one  after 
another  (as  seen  in  Fig.  3),  considering  their  constraints  and  lim¬ 
itations  on  the  storable  power.  If  any  constraint  is  violated,  P,(t) 
of  the  device  should  be  modified  so  that  the  violated  constraint  is 
satisfied.  These  constraints  violations  or  incompatibilities  and  the 
necessary  modification  to  P,(t)  in  each  case  are,  as  follows: 

Case  (1-1 ):  The  power  assigned  to  device  i,  P,(t),  exceeds  its  rated 
power.  In  this  case,  it  is  set  to  the  rated  power  ( Prated,i )  of  the  device, 
i.e.: 


P?{t)  =  Prated , 


(19) 


Case  (1-2):  P,(t)  exceeds  the  remaining  storage  capacity  of  device  i. 
In  this  case,  the  power  is  injected  to  the  device  until  the  maximum 
allowable  level  of  the  stored  energy  is  reached  (SOCmax): 


rated  J 


(SOCmaxi  -  SOCft  —  1 )) 


(20) 


Case  (1-4):  The  operating  limitations  of  the  device  considering  its 
dynamics  are  not  satisfied.  For  instance,  the  change  in  the  power 
share  of  the  device  from  time  step  t  to  t+ 1  is  more  than  the  maxi¬ 
mum  tolerable  rate  of  the  power  change  due  to  the  system  inertia 
(i.e.,  the  power  ramp  rate,  represented  by  r).  This  inertia  implies  the 
impossibility  of  timely  tracking  very  rapid  input  power  fluctuations 
[23]. 

P*(t)  =  r  •  At  (22) 

In  above  mentioned  cases,  the  power  assigned  to  other  devices, 
Pj(t),  should  also  be  modified  because  the  constraint  did  not  allow 
the  device  i  to  fulfill  the  storage  of  its  initially  assigned  P,(t),  calcu¬ 
lated  in  (18).  Therefore,  a  portion  of  P,(t)  is  still  remained,  namely 
Prem.i(t ),  which  other  devices  may  be  able  to  absorb  it.  The  modifi¬ 
cation  procedure  is  named  Charge  Share  Correction  (CSC),  written 
as  follows: 

Prem,i(t)  =  Pi(t)-P*(t)  (23) 

CS*(t)  =  0  (24) 

CSUt)  =  CSj  (  1  +  CS| -  )  for;  +  i,  j  =  1 , . . . ,  n  (25) 

V  l^k=\,k  +  fSk  J 

i-1 

P*(t)  =  Pj(t )  +  CSJ(t)  '  ^  'Prem ,  k(0  (26) 

k= 1 

where  CS,  is  the  charge  share  of  all  devices  except  device  i,  and  PT  (t) 

J  J 

is  the  modified  power  assumed  for  the  storage  device  j.  Eq.  (24) 
means  that  the  device  i  has  absorbed  all  possible  amount  of  power 
it  could,  and  there  is  no  more  power  stored  in  it  during  the  current 
time  step  t.  If  one  of  the  constraints  corresponding  to  this  device  is 
violated,  the  same  modification  procedure  is  then  performed,  for 
that  constraint. 

For  Prem(t)  <  0:  In  this  case,  the  load  demand  is  partially  supplied 
by  RES.  Hence,  Discharge  Shares  (DS)  should  be  determined  and 
assigned  to  the  back-up  devices.  Utilizing  these  back-up  devices, 
the  effort  is  made  to  fully  compensate  the  remaining  power  by  the 
power  from  suppliers  other  than  the  RES  ( Pbackup  )>  and  try  to  achieve 
P backup  ~  \Prem  I- 

The  power  to  be  drawn  from  backup  device  i,  P,(t),  is  written  as 
follows: 


P«(0  =  DSfk)  ■  Prem(t)  i  =  1 ,  . . . ,  m  (27) 

where  m  is  the  total  number  of  back-up  devices  capable  of  feeding 
power  to  the  load  and  DSfk)  is  the  discharge  share  of  the  back-up 
device  i  in  month  k  (/<=  1 . 12).  All  DSs  fall  within  the  range  (0,1). 

As  shown  in  Fig.  3,  the  status  of  all  back-up  devices  is  checked 
one  after  another  considering  their  constraints  and  limitations  on 
the  power  that  they  are  able  to  serve  in  the  current  time  step.  If 
any  constraint  is  violated,  P,(t)  of  the  device  should  be  modified,  as 
follows,  so  that  the  violated  constraint  is  satisfied. 

Case  (2-1):  P,(t)  exceeds  the  rated  power  of  device  i.  (Eq.  (21) 
should  be  applied.) 

Case  (2-2):  The  stored  energy  in  the  device  in  the  current  sim¬ 
ulation  time  step  is  not  as  adequate  as  to  supply  P,(t).  In  this  case, 
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the  power  is  drawn  from  the  device  until  the  minimum  allowable 
level  of  the  stored  energy  is  reached  (50Cmin  i  ): 

P*(t)  =  %gi(SOC,(t  -  1 )  -  SOCminJ)  (28) 

The  value  of  SOCmin  is  obtained  within  the  context  of  the  opti¬ 
mization  procedure. 

Cases  (2-3)  and  (2-4):  These  cases  are  the  same  as  cases  (1-3) 
and  (1-4),  respectively. 

In  the  above  mentioned  cases,  the  power  assigned  to  other 
devices,  Pj(t),  should  also  be  modified  because  the  constraint  did 


not  allow  device  i  to  fulfill  the  supply  of  its  initially  assigned  P,-(t), 
calculated  in  (27),  and  a  portion  of  P,-(t)  is  still  remained,  namely 
Prem,M,  which  other  devices  may  be  able  to  supply  it.  The  modifica¬ 
tion  procedure  is  named  Discharge  Share  Correction  (DSC),  written 
the  same  as  CSC  (Eqs.  (23)-(25)),  with  only  replacing  all  CS  variables 
by  DS,  and  n  by  m. 

While  any  of  the  above  situations  take  place  and  CSC  or  DSC  is 
performed,  it  is  possible  that  again,  the  modified  dispatch  shares  are 
incompatible  with  system  operation  constraints.  In  this  case,  the 
CSC  or  DSC  is  performed  for  the  second  time  so  that  the  constraints 
are  satisfied.  This  is  the  last  effort  for  the  determination  of  dispatch 
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Fig.  4.  Membership  functions  of  objective  functions  and  constraints. 


shares.  Then  the  remaining  excess  or  deficient  power  is  considered 
as  discarded  or  unmet  power  in  the  current  simulation  time  step, 
respectively  (as  presented  in  Fig.  3). 

4.  Optimization  algorithm 

The  multi-objective  optimization  algorithm  proposed  in  this 
paper,  takes  the  advantage  of  the  differential  evolution  algorithm 
accompanied  with  fuzzy-multi  objective  technique. 

4.2.  Multi- objective  optimization  using  fuzzy  technique 

The  multi-objective  problem  is  generally  solved  by  three  types 
of  methods.  The  first  one  is  pareto-based  approach  to  get  a  set 
of  non-dominated  solutions  in  the  process  of  optimization  [24]. 
The  second  one  is  the  coefficient  method  and  the  last  method  is 
transforming  the  multiple  objective  function  into  a  single  objec¬ 
tive  model  and  optimizing  it  through  single  objective  strategies 
[24].  In  this  method,  the  determination  of  appropriate  values  for 
the  coefficients  is  of  very  important.  In  addition,  the  results  are 
highly  dependent  on  and  sensitive  to  selected  values  for  coeffi¬ 
cients.  In  the  method  developed  by  Bellman  and  Zadeh  [25],  the 
single  objective  problem  is  achieved  by  maximizing  the  minimum 
degree  of  satisfaction  among  membership  functions. 

The  basic  idea  in  fuzzy  optimization  is  to  optimize  objective 
function  and  constraints  simultaneously  [26].  The  fuzzy  decision 
is  marked  out  due  to  the  intersection  of  fuzzy  objectives  and  fuzzy 
constraints.  The  first  operation  is  the  fuzzification  process  of  the 
merged  objective  function  and  the  constraints.  In  this  procedure, 
two  types  of  function  /x(x)  are  defined  for  each  objective  function 
and  constraint,  as  shown  in  Fig.  4.  In  this  figure,  the  minimum  value 
for  each  objective  is  obtained  by  the  single  objective  optimization 
and  the  maximum  value  is  specified  by  the  initial  set  of  control 
variables  in  the  optimization  algorithm.  According  to  this  method, 
it  is  possible  to  change  the  effectiveness  of  any  objective  function 
by  reducing  or  increasing  its  specified  maximum  value.  In  other 
words,  when  the  specified  maximum  value  of  an  objective  function 
decreases,  it  is  considered  to  be  more  important  in  the  optimiza¬ 
tion  than  the  previous  state  and  vice  versa.  These  membership 


functions  are  initially  combined  by  “and”  operator  (minimum).  This 
procedure  can  be  expressed  by  the  following  equation: 

/idM  =  min(/xj1(x),  /x/2(x), . . . ,  /xci(x),  /xc2(x), . . .)  (29) 

where  /xD(x)  represents  the  membership  function  of  the  optimal 
decision  function  and  is  treated  as  the  evaluation  value  in  the  opti¬ 
mization  algorithm. 

The  membership  values  express  the  degree  of  satisfaction  for 
each  objective.  Highly  satisfied  objective  is  given  a  low  value, 
though  lowly  satisfied  one  is  assigned  a  high  value.  Hence,  the 
multi-objective  problem  can  be  transformed  into  the  following 
maximization  problem  subject  to  a  crisp  constraint  set: 

max  //D(x,  u)s.t.  H(x,  u)  =  0,  C(x,  u)  <  0  (30) 

where  /xD(x,u)  is  called  the  fuzzy  index. 

4.2.  Differential  evolution  algorithm 

Differential  Evolution  Algorithm  (DEA),  introduced  by  Storn 
and  Price  [27],  is  a  simple  population  based  stochastic  parallel 
search  evolution  algorithm  for  global  optimization  and  is  capable 
of  handling  non-differentiable,  nonlinear  and  multi-modal  objec¬ 
tive  functions  [28].  In  DEA,  the  population  consists  of  real-valued 
vectors  with  dimension  D,  which  is  equal  to  the  number  of  control 
variables.  The  initial  population  with  the  size  Np,  is  uniformly  dis¬ 
tributed  in  the  search  space,  falling  within  variables’  boundaries. 
The  procedure  of  this  algorithm  is  shown  in  Fig.  5.  The  algorithm  is 
described  in  the  following  steps: 

Step  (1),  Initialization 

For  each  control  variable  k  with  lower  bound  xj^in  and  upper 

bound  xfnax,  initial  values  are  randomly  and  uniformly  chosen  in 
the  interval 

XtG  =  Xmin  +  rand  [°’1]><(xmax-Xmin)  *  e  [1,  Wp],  fce[l,D](31) 

Step  (2),  Evaluation 

In  DEA,  after  generation  of  the  initial  population  XG,  all  of  its 
consisting  vectors  are  evaluated  by  the  calculation  of  the  objective 
function  in  its  subprogram,  and  then  the  vectors  are  sorted  accord¬ 
ing  to  their  objective  function  value.  In  the  present  work,  DEA  is 
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Modification  &  constraints  verification  (step  3-3) 
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Replace  the  current  particle  with  the  trial  vector 
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Fig.  5.  Procedure  of  the  used  DEA. 


The  elements  of  the  mutant  vector  enter  the  trial  vector  with  prob¬ 
ability  CR  as  follows: 

,,  ,  f  Vj.tc+1  if  rand,, ,  <  CR  or j  =  /ran  d  f331 

j,l,G+l  |  x.  .  c  jf  randj  .  >  CRandj  ^  /fand 

where  randy,  e  [0,1],  /rand  is  an  integer  randomly  selected  from 
[1,2,. . .,  D].  7rand  ensures  that  one  of  the  trial  vectors  is  selected 
among  the  mutant  vectors. 

Step  (3-3),  Modification  and  constraints  verification : 

This  step  prepares  the  trial  vectors  for  the  next  step  by  veri¬ 
fying  the  variables  whether  they  meet  their  constraints,  given  in 
Sections  3.2  and  4.4.  Otherwise,  for  the  vectors  containing  vari¬ 
ables  with  unsatisfied  constraints,  mutation  and  crossover  steps 
are  repeated  for  a  specified  number  of  times  or  until  the  problem 
is  fixed.  Nevertheless,  these  vectors  will  be  regenerated,  similar  to 
step  (1). 

Step  (3-4),  Selection : 

The  vectors  generated  in  the  previous  step  are  evaluated  and 
sorted  as  described  in  the  step  (2).  Then,  the  trial  vector  U,G+1  is 
compared  with  the  corresponding  vector  Xi  G  and  the  one  with  the 
better  fitness  value  is  admitted  to  the  next  iteration  as  XliG+1 : 


*;,c+ 1 


(  uuc+ 1 

1 


if/awi)  <m,c)  i€[1  N  ] 

otherwise  ’  p 


(34) 


43.  Representation  of  uncertainty 

Natural  characteristics  of  wind  and  solar  energy  impose  uncer¬ 
tainty  in  their  design  and  operation.  Hence,  considering  various 
possible  scenarios  in  the  model  of  these  resources  can  lead  to 
more  realistic  decisions.  The  uncertain  parameters  are  expressed 
by  probability  distributions,  showing  the  range  of  values  that  a 
parameter  could  take,  and  also  accounting  for  the  probability  of  the 
occurrence  of  each  value  in  the  considered  range  [12].  In  this  paper, 
Weibull  and  Beta  PDFs  are  used  for  the  wind  speed  and  the  solar 
irradiation,  respectively  [21].  Using  this  solution,  seven  scenarios 
are  generated,  to  model  the  uncertainty  of  these  resources. 


accompanied  with  fuzzy  technique  to  handle  the  multi-objective 
problem.  Hence,  this  step  has  some  differences,  described  in  Section 
4.4. 

Step  (3),  Iteration 

This  step  is  the  main  procedure  of  DEA  to  conduct  the  final 
results,  which  is  performed  iteratively,  until  either  an  acceptable 
solution  has  been  reached,  or  for  at  least  a  specified  number  of  iter¬ 
ations,  no  further  improvement  in  the  solution  is  obtained,  or  a 
predefined  number  of  iterations  are  completed  [28]. 

Step  (3-1),  Mutation : 

In  DEA,  the  perturbations  applied  to  the  current  population  to 
generate  new  population  are  derived  from  the  current  population 
itself,  and  no  predefined  probability  density  function  is  considered. 
For  a  given  vector  of  control  variables,  XI>G,  in  the  population,  two 
vectors  (Xr2  G  andXrs?G)  are  randomly  selected  such  that  the  indices 
i,  r2  and  r3  are  mutually  different.  Then,  the  weighted  difference 
between  two  vectors  is  added  to  the  vector  with  the  best  fitness 
function,  Xb  G,  written  as  follows: 

^i,G+l  =^b,G •  (^r2,G -Xr3?G)  r2 ,  r3  G  {1 , 2,  .  .  .  ,  NP]  (32) 

where  F,  a  constant  from  (0,  2),  characterizes  the  amount  of  the 
movement  of  vectors  within  the  search  space. 

Step  (3-2),  Crossover : 

In  this  step,  the  trial  vector  Ul>G+1  is  developed  from  the  elements 
of  the  target  vector  X,G  and  the  elements  of  the  mutant  vector  VI>G+1 . 


4.4.  Implementation 

A  program  in  MATLAB  environment  is  developed  to  determine 
the  best  set  of  design  and  PMS  parameters.  Furthermore,  the  con¬ 
vergence  and  computational  time  for  the  algorithm  is  evaluated. 
The  same  system  is  simulated  via  two  other  optimization  methods 
widely  applied  to  renewable  energy  [29],  namely  Genetic  Algo¬ 
rithm  (GA)  [30]  and  Linearly  Decreasing  Inertia  Particle  Swarm 
Optimization  (LDI-PSO)  [31,32]  and  their  corresponding  results  are 
compared  in  the  Section  5.2.  The  performance  of  the  proposed  PMS 
is  also  evaluated  by  comparison  with  another  PMS. 

There  are  two  main  modules  in  the  developed  program,  the 
optimization  module  and  the  system  simulation  module.  The  first 
module  generates  and  modifies  a  valid  vector  of  control  variables 
to  be  evaluated  by  the  system  simulation  module  in  each  iteration. 
When  the  latter  module  is  invoked,  the  input  vector  is  processed. 
The  module  simulates  the  operation  of  the  system  for  each  scenario, 
calculates  the  operating  point  of  each  system  component  and  deter¬ 
mines  ptffix)  and  /zCI(x),  and  the  resultant  /xD(x),  using  the  resultant 
values  for  objective  functions.  These  values  are  transformed  into  a 
single  value  using  the  described  fuzzy  technique  and  become  the 
evaluation  value  for  the  considered  scenario.  The  final  evaluation 
value  for  an  input  vector  is  the  average  of  values  obtained  by  using 
data  of  all  generated  scenarios  [12],  which  is  then  admitted  back  to 
the  optimization  module. 

Therefore,  a  valid  control  vector  is  obtained  by  assigning  val¬ 
ues  to  the  control  variables  that  lie  within  allowed  bounds  of  each 
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individual.  However,  some  modifications  to  the  original  optimiza¬ 
tion  algorithm  are  made  to  match  the  specific  characteristics  of  the 
present  problem,  as  follows: 

-  Some  of  the  control  variables,  including  number  and  type  of  each 
component,  must  be  non-negative  integers.  Therefore,  the  opti¬ 
mization  procedure  should  handle  the  discrete  variables  in  all 
initialization,  crossover  and  mutation  steps. 

-  The  upper  limit,  associated  with  CSj  and  DSj  (i  =  2,. .  .,n)  float  in 
accordance  with  the  values  of  CSk  and  DSk  (/<  =  1,. .  .,i  -  1 ),  because 
CSj  and  DSj  must  satisfy  the  equality  constraint  Ylk=t = 

Y^k=\ =  1  •  This  constraint  modification  is  applied  to  the  value 
that  the  optimization  procedure  assign  to  CSj  and  DSj,  from  device 
i  =  1  through  n,  one  by  one,  i.e.: 

i—i 

CS™*  =  1 ,  CS™*  =  1  -  y)csk  (i  =  2 . n)  (35) 

k= l 

A  similar  equation  is  used  by  replacing  CS  by  DS  and  n  by  m. 

-  The  integration  of  the  fuzzy  model  with  the  optimization  pro¬ 
cedure  for  handling  the  multi-objective  problem  consists  of 
developing  membership  functions  of  objectives  as  well  as  con¬ 
straints  imposed  on  the  operation  of  the  system.  By  solving 
a  single  objective  optimization  for  each  objective  function,  its 
minimum  and  maximum  values  are  obtained  and  the  corre¬ 
sponding  membership  function  is  achieved.  By  the  way,  using 
maximum  and  minimum  power  ramp  rates  of  each  compo¬ 
nent  ( ramp  rateimin  and  ramp  ratei<max),  rated  limitations  on  SOQ 
(50Cj  min  and  SOQ  max)  and  power  boundaries  of  each  component 
( rated  powerim[n  and  rated  poweq  max),  all  as  the  operation  con¬ 
straints,  their  corresponding  membership  functions  are  obtained 
(as  shown  in  Fig.  4). 

-To  treat  the  uncertain  parameters,  several  scenarios  are  generated 
and  simulated.  Hence,  all  the  three  objective  functions  introduced 
in  Section  3.1  are  calculated  for  all  the  scenarios.  The  algorithm 
utilizes  the  concept  of  the  “expected  objective  function”,  which 
aggregates  all  the  objective  function  values  corresponding  to 
the  simulated  scenarios  representing  the  uncertain  parameters, 
in  the  form  of  an  average  value  [12].  Then,  using  the  expected 
objective  function  values,  the  fuzzy  index  is  calculated,  and  the 
algorithm  continues  with  the  same  algorithmic  operations  found 
in  DEA.  In  this  way,  the  results  represent  the  solution  with  an 
overall  consideration  of  all  scenarios  associated  with  the  uncer¬ 
tain  parameters  in  the  system. 


5.  Case  study  and  results 

5.1.  System  data 

The  hybrid  energy  system  under  study  has  been  described  in 
Section  2.  In  practical  system  planning,  the  designer  has  to  choose 
the  system  components  from  a  set  of  commercially  available  ones. 
As  a  case  study,  various  models  and  types  of  components  with  dif¬ 
ferent  rated  powers,  costs  and  other  specification  [  1 5 ],  are  the  input 
to  the  design  procedure.  The  system  location  is  a  region  in  Arde- 
bil  city  (38°14/57//N/48°18/5//’E),  in  the  north  west  of  Iran.  For  this 
region,  the  hourly  data  for  solar  irradiance,  wind  speed  and  ambient 
temperature  all  provided  within  one  year  period  [2].  The  consid¬ 
ered  load  profile  is  the  IEEE  reliability  test  system  [33]  with  30  kW 
peak  power. 

5.2.  Results 

The  data  are  analyzed  to  determine  the  optimum  design  and 
PMS  parameters  for  the  case  study.  The  results  of  the  used  opti¬ 
mization  algorithms,  i.e.,  GA,  LDI-PSO  and  DEA,  are  presented  in 
Table  1.  It  can  be  seen  in  this  table  that  DEA  revealed  the  best  per¬ 
formance  with  respect  to  optimization  of  objective  functions  and 
the  resultant  fuzzy  index.  Hence,  DEA  is  nominated  to  demonstrate 
the  results  of  this  problem  in  more  detail.  Since,  this  problem  should 
be  solved  in  off-line  mode,  the  computation  time  and  convergence 
speed  are  not  of  great  concern. 

In  order  to  well  present  the  applicability  and  efficacy  of  the  pro¬ 
posed  PMS,  which  is  the  significant  idea  of  this  paper,  a  simulation 
is  also  carried  out  by  using  another  PMS,  referred  to  as  PMS  (b). 
Based  on  the  PMS  (b)  which  is  analogous  with  the  PMS  (dispatch 
strategy)  employed  in  the  HOMER  software  [8],  Prem(t)  in  time  step 
t  is  dispatched  among  components  by  the  priorities  defined  by 
designer,  unless  the  marginal  cost  of  energy  storage  or  production 
of  the  dispatched  power  by  the  components  is  in  contrast  with 


Table  1 

Results  of  three  evolutionary  algorithms. 


PMS  (a) 

GA 

LDI-PSO 

DEA 

PMS  (b) 

DEA 

Cost  (€) 

1,042,242 

967,375 

912,560 

918,337 

Fuel  emissions  (kg) 

2086.2 

2632.8 

1827.7 

2888.6 

LPSP  (W/h) 

0.0067 

0.0083 

0.0072 

0.0075 

Number  of  iterations  to  converge 

29550 

196 

227 

153 

Computation  time  (s) 

173,520 

121,311 

134,548 

72,549 

Fuzzy  index 

0.9980 

0.9973 

0.9996 

0.9915 
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Table  2 

Optimized  values  for  design  parameters  of  the  case  study. 


Design  data  DEA 


Equipments  data 

PV  panel 

Wind  turbine 
Battery 

Fuel  cell 

Electrolyzer 

Diesel  generator 
Hydrogen  tank 

Inverter  nominal  power  (kVA) 

Installation  data 

Wind  tower  height  (m) 

PV  panels  monthly 
tilt  angle  -$u  (°) 


Number  9 

Type  10 

Number  6 

Type  6 

Number  28 

Type  5 

Number  6 

Type  4 

Number  8 

Type  4 

Number  7 

Type  3 

Number  6 

Type  35 


32 

67 

58 
63 
46 
38 
31 
22 
30 
26 
43 

59 
53 


the  predefined  priorities.  In  this  case,  the  priority  of  components 
to  store  or  produce  Prem(t)  is  sorted  ascendingly  according  to 
their  corresponding  cost  of  energy  storage  or  production.  The 
optimization  results  using  PMS  (b)  are  included  in  Table  1,  too.  It 
can  be  observed  that  PMS  (a)  have  better  results  than  PMS  (b).  The 
fuel  emission  is  comparably  less  than  that  of  PMS  (b),  showing 
that  the  utilization  of  back-up  and  storage  devices  other  than  the 
diesel  generator  is  more  effective  when  PMS  (a)  is  employed.  Fig.  6 
presents  the  LOEE  in  both  PMS  cases,  confirming  the  better  LPSP  of 
Table  1  for  the  case  of  PMS  (a). 

Table  2  presents  the  result  of  design  parameters,  including  the 
type  and  number  of  each  component.  The  installation  data  con¬ 
sisting  of  the  monthly  PV  panels  tilt  angle  and  wind  turbines 
installation  height  are  also  listed  in  the  table. 

The  monthly  values  of  charge  and  discharge  shares  are  pre¬ 
sented  in  Figs.  7  and  8,  respectively.  Fig.  7  presents  that  the  DS 
of  fuel  cells  is  more  significant  in  winter’s  months  than  summer’s 
months  and  the  opposite  is  true  for  diesel  generators.  The  values 
obtained  for  DS  parameters  also  have  great  influence  on  the  CS 


4 


7 


— I—  Battery 
Diesel 
Fuel  cell 


Fig.  7.  Monthly  dispatch  shares  of  backup  devices  for  providing  the  deficient 
demand. 


- Battery 

X  1  Electrolyzer 


Fig.  8.  Monthly  dispatch  shares  storage  devices  for  saving  the  excess  power. 

Table  3 

Optimized  values  of  limitations  of  storage  devices. 


Storage  device  limitations  DEA 


Minimum  state  of  charge  of  batteries  (SOCmin  %)  21 

Maximum  state  of  charge  of  batteries  (SOCmax  %)  100 

Minimum  level  of  Hydrogen  in  tanks  (%)  33 


values.  In  fact,  CS  and  DS  parameters  are  interdependent  for  compo¬ 
nents  that  utilize  the  same  energy  storage  device  in  both  charging 
and  discharging  states.  For  instance,  electrolyzers  and  fuel  cells  use 
the  same  energy  container  namely  hydrogen  tanks  to  store  energy 
and  consume  it,  respectively.  Table  3  contains  optimized  values  of 
operating  limitations  of  storage  devices. 

The  attained  values  for  parameters  of  PMS  by  represent  nearly 
the  overall  optimum  operating  condition  of  the  system  considering 
a  compromise  among  all  the  considered  objectives.  The  parameters 
are  separately  calculated  for  each  month,  representing  the  influ¬ 
ence  of  monthly  and  seasonal  changes  in  load  demand  and  climatic 
patterns  on  the  utilization  of  the  hybrid  system.  Hence,  the  results 
are  adapted  with  any  circumstances  that  alter  the  operation  of  the 
system  resulting  in  more  accuracy  and  optimality. 

6.  Conclusions 

This  study  has  been  dedicated  to  review  optimal  design  and 
operation  strategy  selection  of  hybrid  RES-based  stand-alone 
energy  systems,  including  various  generators  and  storage  media. 
Considering  resource  uncertainties,  associated  with  wind  speed 
and  solar  irradiation,  a  novel  method  has  been  proposed,  to  deter¬ 
mine  the  power  management  strategy  of  the  system  along  with 
sizing  parameters  of  system  components  so  that  the  overall  cost 
of  the  system  during  its  useful  lifetime,  unmet  load  and  pollutant 
emissions  have  been  minimized.  The  PMS  parameters  are  monthly 
charge  shares  (CS)  of  storage  devices  and  monthly  discharge  shares 
(DS)  of  generator  devices,  limitations  on  their  start-up/shut-down 
power  thresholds  (HBR)  and  on  their  level  of  available  energy  (SOC). 

Based  on  the  proposed  method,  the  values  of  PMS  parameters, 
the  sizing  and  design  parameters  have  been  determined  by  the 
iterative  optimization  algorithm.  The  PMS  parameters  have  been 
separately  determined  for  each  month  to  adapt  monthly  variations 
in  the  load  and  climate  patterns.  Then  in  each  iteration,  these  values 
have  been  applied  to  the  hourly  simulation  of  the  system  opera¬ 
tion  and  evaluated  in  each  single  time  step  to  meet  the  operational 
constraints  of  the  components,  consisting  of  the  nominal  power, 
SOCmin  and  SOCmax,  HBR,  and  power  ramp  rate  of  each  device.  Oth¬ 
erwise,  these  values  have  been  modified  by  either  of  two  operators, 
namely  CSC  or  DSC,  in  that  time  step.  These  modifications  have 
ensured  that  while  the  optimum  operation  strategy  and  commit¬ 
ment  of  components  has  been  attained,  all  constraints  associated 
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with  operational  characteristics  of  the  components  have  been  sat¬ 
isfied.  However,  if  after  completion  of  all  steps  of  simulation  in  a 
time  step,  the  constraints  were  still  unsatisfied,  some  discarded  or 
unsupplied  energy  remains.  The  summation  of  these  values  during 
all  time  steps  of  the  simulation  has  been  involved  in  the  objective 
functions  to  be  minimized. 

The  proposed  method  has  tested  on  a  hybrid 
PV-wind-diesel-hydrogen-battery  system.  To  efficiently  handle 
the  nonlinear  mixed-integer  multi-objective  optimization  prob¬ 
lem,  a  modified  version  of  DEA  accompanied  by  fuzzy  technique 
has  implemented.  The  system  has  been  optimized  for  summations 
of  all  objective  values  of  each  function  calculated  for  all  scenarios 
generated  by  applied  uncertainty  models.  The  presented  results  of 
the  solution  can  be  categorized  into  three  groups: 

-  Size  (design)  data:  The  optimum  number  and  type  of  each  com¬ 
ponent  of  the  system  to  be  installed. 

-  Installation  data:  The  optimum  values  for  monthly  inclination 
angles  of  solar  panels  and  the  optimum  tower  height  for  wind 
turbines  installation  to  maximize  the  exploitable  power  from  RES. 

-  Operation  strategy  (PMS):  The  monthly  CS  and  DS,  and  boundaries 
associated  with  SOC. 

By  comparing  the  results  of  design  parameters  and  objective 
functions  values  (accompanied  by  that  of  the  system  operation 
and  unsupplied  energy),  with  the  results  of  applying  a  user- 
defined  PMS,  it  has  been  concluded  that  the  design  decisions  and 
operation  of  such  systems  are  highly  affected  by  the  employed 
PMS;  and  the  superiority  of  the  proposed  PMS  has  been  con¬ 
firmed.  It  should  be  mentioned  that  the  proposed  method  can 
be  effectively  employed  for  any  composition  of  hybrid  energy 
systems. 
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